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ABSTRACT 

In order to study the outflows from accretion disks, we solve the set of hydrody- 
namic equations for accretion disks in the spherical coordinates (rOcp) to obtain the 
explicit structure along the 6 direction. Using self-similar assumptions in the radial 
direction, we change the equations to a set of ordinary differential equations (ODEs) 
about the ^-coordinate, which are then solved with symmetrical boundary conditions 
in the equatorial plane, and the velocity field is obtained. The a viscosity prescrip- 
tion is applied and an advective factor / is used to simplify the energy equation. The 
results display thinner, quasi-Keplerian disks for Shakura-Sunyaev Disks (SSDs) and 
thicker, sub-Keplerian disks for Advection Dominated Accretion Flows (ADAFs) and 
slim disks, which are consistent with previous popular analytical models. However, an 
inflow region and an outflow region always exist, except when the viscosity parameter 
a is too large, which supports the results of some recent numerical simulation works. 
Our results indicate that the outflows should be common in various accretion disks and 
may be stronger in slim disks, where both advection and radiation pressure are dom- 
inant. We also present the structure dependence on the input parameters and discuss 
their physical meanings. The caveats of this work and possible improvements in the 
future are discussed. 

Subject headings: accretion, accretion disks - hydrodynamics - black hole physics 



1. INTRODUCTION 



The accretion disk models have been developed much in the past several decades. Many disk 
models have been proposed and some of them are widely adopted in astrophysical studies nowa- 
days, including but not restricted to Shakura-Sunyaev Disk (SSD; see Shakura & Sunyaev 1973), 



'E-mail: chengliang.jiao® pku .edu.cn 
-E-mail: wuxb@pku.edu.cn 



-2- 



Advection Dominated Accretion Flow (ADAF; see Narayan & Yi 1994; Abramowicz et al. 1995), 
and slim disk (Abramowicz et al. 1988). However, the outflow structure of accretion disks still 
remains an open problem. The equations that describe the hydrodynamic processes are the Navier- 
Stokes equations, which are quite difficult to solve in the case of accretion disks which involve 
viscosity and radiation. Therefore, in most works, some kind of simplification, such as one-zone 
or polytropic distribution and hydrostatic equilibrium, are usually applied in the vertical direction, 
and the vertical variation (z dependence in cylindrical coordinates) of the velocity field is usually 
neglected. In this way the equations are changed to a set of ordinary difi'erential equations (ODEs) 
in the radial direction, which can be solved numerically. However, by taking these assumptions, 
one cannot get a clear picture of the vertical structure of accretion flows; the velocity is always 
radially inward and no mass will cross the disk surface, displaying no outflow structure. Among 
the exceptions is a work done by Narayan & Yi in 1995(hereafter NY95), which used self-similar 
assumptions in the radial direction and solved the structure along the 6 direction in spherical coor- 
dinates (rOc/)). However, in their work they assumed vgi = and thus cannot get a proper velocity 
field, and their solutions compose of only pure inflow. They argued that the Bernoulli parameter 
is positive in their solutions so that a bipolar outflow is expected to develop near the vertical axis. 
Blandford & Begelman (1999, hereafter BB99) relaxed the mass conservation assumption and 
assumed that the mass inflow rate varies with radius, and obtained solutions with outflow called 
adiabatic inflow-outflow solutions (ADIOS). Their solutions are one-dimensional self-similar solu- 
tions that are height- averaged and they also applied the Bernoulli parameter to argue the presence 
of outflow. However, Abramowicz et al.(2000) pointed out that the positive Bernoulli function 
(instead of Bernoulli parameter, which is defined in convenience of self-similar assumptions) is 
not sufficient for outflow (also see the simulation works done by Stone et al. 1999 and Yuan & 
Bu 2010). Blandford & Begelman (2004) furthered their work and presented some self-similar 
two-dimensional solutions of radiatively inefficient accretion flows with outflow. They assumed 
hydrostatic equilibrium in the vertical direction and that convection dominates the heat transport, 
which may only be applicable in certain cases. Xu & Chen (1997, hereafter XC97) relaxed = 
and obtained two types of solutions with outflow: accretion and ejection solutions. However their 
solutions require the net accretion rate to be 0, which is not realistic. Xue & Wang(2005, hereafter 
XW05) followed NY95 and solved the disk structure along the 9 direction considering vg. They 
arbitrarily set a disk surface, at which = 0, and the sound speed on the surface for their calcula- 
tion. Their solutions display a field of inflow near the equatorial plane with wind blowing out of 
the upper boundary, however the boundary is set as an input parameter rather than being calculated, 
and they only investigated several cases of ADAFs. In §3.2 we will see that according to our calcu- 
lation, their assumption doesn't hold for accretion flows with large a value. Sadowski et al.(2010) 
abandoned the self-similar assumptions and solved the accretion disk structure in the radial and 
vertical direction simultaneously. Because the Navier-Stokes equations for accretion disks cannot 
be decoupled intrinsically, they adopted other assumptions, e.g. the disk is not geometrically thick. 
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to derive the equations. As in their work they didn't consider while Vr and were supposed 
not to vary vertically, they were not able to study outflows. In summary, in the analytical model 
study, the vertical or ^-direction structure of accretion disks and outflows can still not be dealt with 
satisfactorily, and many improvements are still needed to be made. 

On the other hand, observationally there are more and more evidences for outflows in accre- 
tion systems, such as Sgr A*(Marrone et al. 2006; Xie & Yuan 2008), soft X-ray transients(Loeb 
et al. 2001) and quasars with blueshifted absorption lines(e.g., PGl 1 15+80; Chartas et al. 2003). 
Many numerical simulation works have also discovered outflows in their results (e.g. Stone et al., 
1999; Igumenshchev & Abramowicz, 2000; Okuda et al. 2005; Ohsuga et al.2005; Ohsuga & 
Mineshige 2007; Ohsuga et al. 2009) . The common existence of outflows in these works inspire 
us to investigate the vertical structure of accretion flows explicitly, to find solutions which can deal 
with vg and positive to get a clear velocity field, and with more reasonable boundary conditions. 
As a first step, we followed the work done by NY95 and XW05, using self-similar assumptions 
in the radial direction and solving the ODEs along the 6 direction in spherical coordinates (r9^). 
We used the a viscosity prescription and assumed that the r0 component of the viscosity stress 
tensor is dominant. By neglecting other components of the viscosity stress tensor, the necessary 
number of boundary conditions is reduced and we will only need the boundary conditions in the 
equatorial plane, which is obatined by symmetry and is thus quite physical. As we didn't set other 
arbitrary restrictions to Vr and vg other than the self-similar assumptions, we can get a velocity field 
containing positive v^, to discuss the flow structure with possible outflows. These assumptions are 
applicable to different kinds of accretion disk models, ranging through SSD, ADAF and slim disk, 
so we also did many calculations with different sets of parameters (which is very difficult to do in 
a numerical simulation work as it would be too time-consuming), trying to find the flow structure 
dependence on the parameters, in order to understand the inflow/outflow mechanism more phys- 
ically. These results can also be helpful to future numerical simulation works when they set the 
input parameters. 

It should be noted that there is another branch of researches which investigate accretion flows 
with an outflow (usually wind) plus accretion disk model. In these researches, the configuration 
of the outflow and accretion disk is usually preset, and the calculations are focused on either the 
outflow or the accretion disk while the influence of the other part is simplified or parameterized 
(e.g. Fukue, 1989; Takahara et al., 1989; Kusunose, 1991; BB99; Misra & Taam, 2001; Fukue, 
2004; Xie & Yuan, 2008). Recently there are also several works done in this way which deal with 
the outflow and the accretion disk simultaneously (Kawabata & Mineshige, 2009; Dotan & Shaviv, 
2010). Compared with our work, in these studies the accretion disk is usually height-integrated and 
the configuration of the accretion flow is assumed rather than being calculated, while in our work 
we solve the full hydrodynamic equations to get the configuration of the accretion flow. Our work 
focuses on studying the general structure of disks, outflows and the physical mechanism behind 
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them, and the results can be complementary to one another. 

In §2 we present the basic equations and assumptions we used in our calculations. In §3.1 we 
discuss our numerical methods and present solutions corresponding to typical parameters of SSD, 
ADAF and slim disk. In §3.2 we show the disk structure dependence on different parameters and 
discuss their physical meanings. In §4 we present our summary and discussion. 



2. EQUATIONS AND ASSUMPTIONS 

We consider the hydrodynamic equations of an accretion flow in the spherical coordinates 
{rO(j)). The flow is assumed to be steady(5/5? = 0) and axisymmetric(5/50 = 0). Thus the conti- 
nuity equation can be written as 

\ d . Id 

^ ^(^pvr) + ;^(sin Opve) = 0. (1) 
H or r sm 6 od 

We assume that for the accretion flow, only the r0-component of the viscous tensor, tr^, is dom- 
inant. And we use the Newtonian gravitational potential, O = -GM/r. Then the equations of 
motion can be written as(Kato et al. 2008, also see Appendix A of XW05) 

dVr Vg dVr . GM I dp 
Vr-^ + —{-TT -Ve) = ^ — (2) 

or r o6 r r^ p or 



dvg 

+ — 

or 



= ^ 

r oU r pr o& 



dr r do r pr^ dr 

In our calculation, we adopt the a prescription of viscosity tr^ = -ap and p is the total pres- 
sure. The shearing box radiation Magnetohydrodynamics (MHD) simulations done by Hirose et 
al. (2009) found that the vertically integrated stress is approximately proportional to the vertically 
averaged total thermal (gas plus radiation) pressure. In our work, however, we also apply this 
relation locally. 

Following NY95, we use the advective factor, / = Qadv/Qvis. i-e. a fraction / of the dissipated 
energy is advected as stored entropy and a fraction (1 -/) is lost due to radiation. In our calculation 
we assume that / is constant in the accretion flow, so the energy equation can be written as 

de vede p dp vedp 5 

P(Vr-^ + -^) - -(Vr-^ + = Pr,l,r—(—), (5) 

dr r d& p dr r da dr r 
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where e is the internal energy of the material and can be expressed as(Kato et al. 2008) 

pe = + 3prad, (6) 

y- 1 

where y is the heat capacity ratio and is considered to be a constant input parameter in our calcu- 
lation. 

We adopt the self-similar assumptions in the radial direction, therefore we seek a solution of 
the form 

P = p(e)r-\ (7) 

Vr = vmJ—, (8) 



ve = Ve(e)J—, (9) 



= v^e)^^, (10) 

p = p{e)GMr"-\ (11) 

This set of self-similar solutions is similar to that of NY95 and is actually identical to that of 
XW05. Compared with those works, here we use the total pressure p in Equation (11) instead of 
the sound speed c,, which is proportional to ^jp|p. As stated in NY95, the only length scale in 
the problem is r and the only frequency is O.^, thus all velocities must scale with radius as rO.^^. In 
NY95, Narayan and Yi set n = 3/2, which implies = according to the continuity equation, thus 
they intrinsically set no outflow for the accretion disk. Here we relax this parameter n following 
BB99 and XW05, to allow outflows from the disk. 

If we substitute Equations (7)-(l 1) into Equations (l)-(4), the r components can be eliminated, 
leaving only the dimensionless fuctions p{6), VriO), Vo(6), v^{0) , p{0), the variable 6 and some 
constants which are set as input parameters. The energy equation (Equation (5)) is somewhat more 
complex. As the internal energy depends diff'erently on gas pressure and radiation pressure, we 
discuss the energy equation more carefully here. First, one may notice that only the total pressure 
p appears in Equations (l)-(4). This is not hard to understand as the dynamical processes don't 
recognize whether the pressure is from gas or radiation, pgas and prad do aff"ect the energy equation 
in a diff'erent way, and to describe this eff'ect, we define the gas pressure ratio yS, 

P^^= . (12) 

P Pgas + Prad 

There is one particular case that when y = 4/3, no matter what the value of is, the solution 
remains the same, as in this case Equation (6) becomes pe = 3p; if the accretion flow is radiation 
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pressure dominated {fi 0), then Equation (6) also becomes pe = 3p, and the flow structure 
won't depend on y; if the accretion flow is gas pressure dominated, then Equation (6) becomes 
pe = pliy - 1), and the result will depend on the value of y. For a general case. Equation (6) can 
be written as 



pe = + 3prad = 

7-1 



+ 3(1 -,8) 



(13) 



.7-1 

Physically ;S is a value between and 1, so pe is always between 3p and pl{y - 1). Then we can 
expect that a general solution should be somewhere between two extreme cases: the gas pressure 
dominated case and the radiation pressure dominated case. For these two extreme cases, we can 
calculate the exact solutions of the problem, as Equation (6) will be simplified to forms without 
p. These two cases also provide the typical solutions we emphasize in §3.1. However we also 
want to know how the solutions change from the gas pressure dominated case to the radiation 
pressure dominated case. For the accretion flows in which both gas pressure and radiation pressure 
are important, we have to assume the gas pressure ratio p to be constant to solve the problem. 
Solving the problem with variant p is beyond the capability of this self-similar treatment, and 
we're planning to do that in our next work. 

As we already mentioned, the radiation pressure dominated case has the same result as that of 
the gas pressure dominated case with y = 4/3, so the question eventually becomes examining how 
the solution changes according to y for the gas pressure dominated case. We define an equivalent 
y here 

Teau = — + 1 = + 1, (14) 

^^"^ pe ^ + 3(l-^)(y-l) 

so that Equation (6) becomes 

pe = — (15) 

Tequ - 1 

and this 7equ represents the equivalent y that the accretion flow would have if the flow is treated as 
gas pressure dominated (even p 1). We can see that if the flow is radiation pressure dominated 
(/S = 0), then no matter what the value of y is, yequ is always 4/3. Here we treat y as a constant, 

and for any specific case, we can get a constant yequ and all the situations can be treated as a gas 
pressure dominated flow with yequ as an input parameter. Usually y is taken as a value between 
7/5(the value for diatomic ideal gas) and 5/3(the value for monatomic ideal gas), so that yequ ranges 
between 4/3 and 5/3. 

With Equations (7)-(l 1) and (15), the Equations (l)-(5) can be reduced to a set of five ODEs: 

= 0, (16) 



dpiO) 
do 



(3 - 2n)vr(e) + 2 1 cot9ve(e) + 
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i{n + \)p{e)+p{e) 



vxef + 2 1 -1 + ve{ef + v^{ef - vM^^ 



= 0, (17) 



-2 cot Ov^ief + ve{e) I vM + j 



0, 



(18) 



- 2{n - 2)ap{e) + p(9) 



Vr(e)v^(9) + 2vg(e) [cot ev^i9) + 



de } 



0, (19) 



2p(e)ve(e)^ + p(e) |p(e) [2(nr^, -n- \)vm - 3af(r^, - l)v^(e)] - 2requv,(e)^| = 0, 

(20) 

with 5 dimensionless fuctions p(6), Vr{9), vg(9), v^(9) , p(9), the variable 9 and four input parameters 
(Q',/,7equ7'^)- This Set of ODEs can be numerically solved with proper boundary conditions. We 
assume the structure of the disk is symmetric to the equatorial plane, and thus we have 

= gO'^:,^ = ± = ± = ^ = ±l=O (21) 
^ d9 d9 d9 d9 ^ ^ 

in which only four conditions are independent. For the last boundary condition we set p(90°) = 1, 
which can be normalized by a scale factor if the effective accretion rate at a certain radius is set 
(NY95, XW05, etc.). These boundary conditions are enough for our calculations and we didn't 
introduce other arbitrary boundary conditions. 



3. NUMERICAL RESULTS 

3.1. Typical Solutions 

We obtained numerical solutions of Equations (16)-(20) with different sets of input parameters 
(a,f,ysqu,n)- Some typical solutions are shown in Figures 1-4. The calculation starts from the 
equatorial plane (9 = 90°) towards the vertical axis (9 = 0°). Generally the mass density p and 
total pressure p will decrease as 9 decreases, and at certain inclination they will get very close to 
as shown in the figures. We take this as the upper boundary of the flow structure. If we continue 
the calculation through this inclination, it will encounter numerical errors. We think the reason is 
that we can't describe the flow near the vertical axis with a simple self-similar solution in the radial 
direction. This eff'ect is also to be expected. If we describe the upper boundary (minimum 9) that 
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we reach in our calculation as 9b, then the effective accretion rate M^s across a sphere at radius r 
within the region calculated by us is 



Meff = 2 I pvr-lnr sin Or ■(7r/lS0°)de = 4nyfGMri-" l v,(%(0)sin0- (;r/180V^, (22) 



which is a function of r unless n = 3/2. In Equation (22) (and subsequently in this paper) negative 
values of M^s represent inflow while positive values represent outflow. If we describe the accretion 
rate in the region between the vertical axis and the inclination as Maxis. then according to the 
steady nature of the flow, we have 



in which M represents the total accretion rate across any sphere at a reasonable radius centered by 
the central accretor and should be a constant for a steady accretion flow. If the solution doesn't 
end at an upper boundary, and instead can describe the entire flow structure in the whole space, 
then ^ = and M^s = M, which should be a constant. According to Equation (22), this can only 
happen in the following two cases: (I) n = 3/2 which enforces M^g not to change with radius 
r; (2) when n 312, the integration term in Equation (22) must be 0, in which case Mgff = 
and is a constant along radius r. The first case was discussed in NY95. Because n = 3/2, r^pvr 
is independent of r, and the continuity equation (1) shows that ve = 0, resulting in a solution in 
which the flow is always radial (with rotation). The second case was discussed in XC97, and the 
fact that M = requires that the outflow rate exactly equals the inflow rate at any radius. However 
this is unrealistic for an accretion flow, which is also discussed in XW05. The reason is that, when 
material is accreted in the form of an accretion flow, gravitational energy is released and part of it 
is changed to internal energy via viscous friction. The restriction that outflow rate equals inflow 
rate requires that the internal energy released from gravitational energy should be fully returned to 
gravitational energy, which violates the second law of thermodynamics. Therefore, the inflow rate 
must be larger than the outflow rate at a certain radius, in order to compensate for the increase in 
entropy in the hydrodynamic process. As we have mentioned before, there are many observations 
of accretion systems with outflows, and they can be neither of the two cases mentioned above (case 
(1) has no outflow while case (2) violates the second law of thermodynamics). So the self-similar 
solutions for an accretion disk with outflow have to be truncated at some inclination. It should 
be noted that this efi'ect has already been discussed in XW05, though they assumed that the self- 
similar solution is only valid for the inflow part. Our solutions, however, show that at least part of 
the outflow (with 6^ < 9 < Oq) can be described with self-similar assumptions. 

From the above analysis, we can also see that n is a very important parameter. When n = 3/2, 
the entire flow structure can be described by the same set of self-similar solutions, however this 
kind of solutions have no outflow. When n 3/2, self-similar solutions can only describe part of 
the entire flow structure. According to Equation (22), when n < 3/2, the eff'ective accretion rate 





Meff + Maxis = M 



(23) 
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Mefr decreases as r decreases, implying that material is lost due to outflows. This kind of solutions 
contain outflows, which is also consistent with the results of many numerical simulation works 
(e.g. Okuda et al. 2005; Ohsuga et al.2005; Ohsuga & Mineshige 2007; Ohsuga et al. 2009). On 
the other hand, when n > 3/2, the solutions have an effective accretion rate M^s that increases 
towards the central accretor, which implies that there must be matter injection into the accretion 
flow from high latitudes. This kind of matter injection is not likely to happen in real cases. So we 
would like to investigate solutions with n < 3/2. Solutions with n = 3/2 have been examined in 
NY95, and they are only for ADAFs, while solutions with n < 3/2 in our work is applicable to a 
much wider range of accretion disk types. So in this subsection, we take n = 1.3 to calculate the 
typical solutions. In this way we can observe how the solutions differ from those of n = 3/2 when 
n doesn't change much. We also calculated how the solutions change with diff'erent n, which will 
be discussed in detail in §3.2. The result shows that, when other parameters are the same, for a 
large range of n with n < 3/2 the solutions display similar structures qualitatively. 

The value of the viscosity parameter a can be inferred both from observations and from MHD 
simulations of magneto-rotational instability (MRI). As reviewed by King et al. (2007), the value 
of a for ionized disks given by numerical simulations is ~ 0.01, while that obtained through 
observations is ~ 0.1 - 0.4. However, the observational determinations of a are strongly model- 
dependent and it still remains an open question what the value of a should be for hot, ionized disks. 
In this subsection we take or = 0.1, which is a typical value used in theoretical models of accretion 
disks. We also discussed how the solutions would change with diff'erent a in §3.2. 

Among the most popular analytical disk models are the SSD, ADAF and slim disk models, so 
here we show the solutions with input parameter corresponding to them respectively. For accretion 
flows composed mostly of fully ionized hydrogen, the material can be regarded as monatomic ideal 
gas, of which y = 5/3, and that is the value we take here (accretion flows composed of unionized 
H2 would have 7 = 7/5. We didn't discuss this situation here, but some relative results can be 
found in §3.2). Figures 1 & 2 are both for the SSD model which has little advection (/ = 0.01), 
while Figure 1 is for the gas pressure dominated region and Figure 2 for the radiation pressure 
dominated region. The ADAF model is treated as advection dominated(/ =1), and gas pressure 
dominated (yS = 1, so that ygqu = 5/3), as shown in Figure 3. The slim disk model is treated as 
advection dominated(/ = 1), and radiation pressure dominated (J3 = 0, so that yequ = 4/3), as 
shown in Figure 4. The other two parameters are all set as or = 0.1 and n = 1.3. In each figure, we 
show the function curves of Vr(9), vg(0), v^{6) , p{6), p{6) and the Mach numbers of Vr and vg, in 
which the adiabatic sound speed is calculated as 



As shown in Equations (7-11), all the velocities scale with ^GMjr. The density p{0) is scaled 




(24) 
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to 1 in the equatorial plane and it can be determined in a real case once the accretion rate at a 
certain radius is given. The pressure p is also scaled, and it satisfies the relation that p(9)/p(6) = 
(p/p)/{GM/r). So once we get the actual density p at a certain radius, we can calculate the cor- 
responding pressure p easily. The Mach number curves for first decreases with 6, reaching at 
some inclination, and then increases. That is because it changes sign along 6 and we use absolute 
values in calculating the Mach numbers. 

There are several common features among Figures 1-4. The value of vq is always non-positive. 
The value of is always negative close to the equatorial plane, and becomes positive at smaller 
inclinations. In each of these figures, there exists a certain inclination at which Vr{9) = 0. This 
^0 has similar meaning with the '^o' in XW05, except that it is obtained through calculation rather 
than being preset as an input parameter in XW05. In the region < < 90°, Vr(0) is negative, 
meaning that the accretion flow is moving towards the central accretor as inflow. So this region 
corresponds to the 'normal' accretion disk part in basic models and we call it the inflow region 
in this paper. In the region 9^ < 9 < ^o, Vr(9) is positive and the accretion flow is moving away 
from the central accretor, and we call this region the outflow region. The region between the upper 
boundary and the vertical axis contains the outflow which blows out of the upper boundary in the 
form of wind as shown in Figure 5, and we call this region the wind region, to distinguish it from 
the outflow region. It should be noted that our solutions actually cannot describe the structure of 
the wind region. However, it is natural to suppose that the flow structure in the wind region with 
< < ^ is also outflow, and the outflow region < ^ < ^o) can either be regarded as a 
transition region between inflow and outflow, or as part of the outflow which can still be described 
with the self-similar assumptions. 

The physical properties shown in Figures 1 & 2 agree quite well with the SSD model. The 
rotation is very close to being Keplerian; the radial velocity is much smaller than the azimuthal 
velocity; and the disk is geometrically thin. However, as mentioned above, our calculation shows 
that even in the SSD case, there exists outflow in the accretion disk. As shown in Figures 5(a) & 
5(b), this outflow is in the form of wind blowing out from the disk surface. Both Vr and vg are 
subsonic in the SSD cases. 

The advection dominated solutions, as shown in Figures 3 & 4, are geometrically much thicker 
than the solutions with low advection, and are able to describe the structure of the majority of 
the space, except only a small region near the vertical axis. The rotation is sub-Keplerian and 
the radial velocity can be comparable with the azimuthal velocity. Both the ADAF and the slim 
disk case show strong outflows compared with the SSD cases, which can be seen more clearly in 
Figure 5. keeps to be subsonic for both cases. However, for the slim disk case (Figure 4), 
becomes supersonic near the vertical axis, displaying a stronger outflow than the other three typical 
solutions. 
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Outflow is the result of the competition among the pressure gradient, the centrifugal force and 
the gravitational force. As we have mentioned before, the key feature of the outflow in our work is 
positive Vr, so we would like to investigate the properties of these forces along the radial direction. 
The radial components of the pressure gradient, the centrifugal force and the gravitational force 
can be written respectively as 

-'-'i = (25) 

p or 

^ = Bief4, (26) 



= 1.2^, (27) 



in which 



AiO) = (28) 
8(9) = vjie). (29) 

We plot the values of A(d), B(9) and A(d)+B(9) for the four typical solutions in Figure 6. The dotted 
lines correspond to the gravitational force which is scaled as 1 ; the dashed lines correspond to the 
radial component of the centrifugal force; the dot-dashed lines correspond to the radial component 
of the pressure gradient; and the solid lines correspond to the sum of the radial components of 
the centrifugal force and the pressure gradient, which drives the outflow. Figure 6(a) shows the 
properties of gas pressure dominated SSDs, while Figure 6(b) shows the properties of radiation 
pressure dominated SSDs. It can be seen that for both of the SSD cases, the influence of pressure 
gradient is very small, and almost in the entire range of the solution, the gravitational force is 
balanced with the sum of the centrifugal force and the pressure gradient in the radial direction. 
This is consistent with the SSD model in which both the advection and the pressure gradient are 
small enough to be neglected. Therefore these two solutions have similar accretion flow structure 
and only some quantitative diff'erences exist between each other. 

On the other hand, for the advection dominated cases, such as the ADAF case (Figure 6(c)) 
and the slim disk case (Figure 6(d)), the pressure gradient plays a more important role than the 
centrifugal force and drives a significant outflow as 6 decreases. The slim disk case has qualitative 
diff'erences from the other three cases. In both SSD cases and the ADAF case, the radial component 
of pressure gradient starts decreasing at smaller inclination angle, and to balance this eff'ect, the 
disk rotates faster as 6 decreases to increase the centrifugal force, as shown in Figures 1-3. However 
in the slim disk case, the radial component of the pressure gradient keeps increasing as 9 decreases. 
It not only drives a much stronger outflow than the other three cases, but also reduces the required 
centripetal force to keep the disk rotating. Thus the disk rotates slower as 9 decreases, as shown in 
Figure 4. 
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As a comparison with NY95, we also present here a solution with parameters a = 0.1, n = 
3/2,yequ = 1.6061and/= linFigureV. This set of parameters corresponds to ' 6/ = 0.rinNY95. 
In this case, the accretion flow doesn't contain outflow and the three components of velocity and 
the sound speed remain constant in the structure. The curves of p{6) and Vr(0) don't agree well 
with NY95, however this is because NY95 also sets the boundary conditions on the vertical axis 
and the values of p(6) and Vr(9) are limited by these boundary conditions. Our solution only uses 
the boundary conditions in the equatorial plane, and can still reach the vertical axis in this case and 
has a reasonable distribution of physical properties, which further proves the applicability of our 
method. 



3.2. Solution Dependence on Input Parameters 

In this section we show how the solutions change according to different input parameters 
respectively, and discuss their physical meanings. As shown in §3.1, the whole space is generally 
divided into three regions: an inflow region, an outflow region and a wind region. As the accretion 
flow is symmetrical both axially and under reflection in the equatorial plane, we only need to 
investigate the flow structure in a quadrant of the r9 plane. In this case, the inflow region starts 
from the equatorial plane (9 = 90°) to the inclination 9 = 9q where Vr{9) = 0; the outflow region 
starts from the inclination 9 = 9Qto the upper boundary of the flow structure at 9 = 9^, and beyond 
that is the wind region (0 < 9 < 9^) in which the motion of the matter doesn't satisfy the self- 
similar assumptions. Here we use 9q and 9^ as indicators of the disk size in 0-direction, which 
directly reflects the size of the inflow region and the whole inflow/outflow region, respectively. In 
the figures showing how the properties of solutions change according to input parameters, each 
line represents a set of solutions, and we assign a unique label for each solution set for clarity. 
Table 1 lists all the solution set labels, their input parameters and corresponding lines at the end of 
§3.2. 

1) Solution dependence on a 

The dependence of 9o and 9^ on a is shown in Figure 8. First we can clearly see that over 
a large range of a, the advection dominated solutions have both larger inflow region and larger 
inflow/outflow region. That's because in these solutions most of the viscous heating is stored in 
the accretion flow as internal energy, resulting in larger pressure at given density and thus larger 
size of the whole inflow/outflow region. 

For solutions with the same advective factor /, gas pressure dominated solutions generally 
have larger inflow region. The inflow region lasts from the equatorial plane to the inclination 9o 
where Vr reaches 0, so its size depends on the behavior of v^. In the equatorial plane, the reflection 
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symmetry gives Equation (21), and when it is substituted into Equations (17), (19) & (20), we can 
get 

piO) o , 
'^{n+\)^+vl{e) + 2{vl{0)-l) = 0, (30) 
p{9) 

-2(n-2)a^ + Vrie)v^ie) = 0, (31) 

2(nrequ - « - l)Vr(e) - 3a/(7equ - l)v^(e) = 0. (32) 

Equation (30) describes the hydrodynamical balance in the radial direction among the pressure 
gradient, the centrifugal force and the gravitational force, which determines the acceleration of Vr 
(the term vjiO) represents the acceleration of Vr under the self-similar assumptions). Equation (31) 
describes how the viscosity affects the angular momentum transfer rate in the equatorial plane. 
Equation (32) describes the energy mechanism of the accretion flow that a fraction / of the dissi- 
pated energy is advected as stored entropy. It should be noted that is generally negative in the 
equatorial plane, so the absolute value of can be written as (according to Equation (32)) 



|v,(90°)| = -v,(90°) = ^v^(90°) 



' -1 



1 - n(yequ - 1) 



(33) 



If v^{90°) remains the same, then in the range of parameters we study in this paper, |vr(90°)| in- 
creases with 7equ- However v^(90°) also changes with input parameters, and the exact value of 
v^(90°) should be calculated via Equations (30)-(32). Figure 10 shows the profile of Vr(90°) vs yequ 
when other parameters are taken as n = 1.3 and or = 0.1. It can be seen that for both the advection 
dominated solutions and the solutions with little advection, |vr(90°)| increases with yequ. With a 
larger value of in the equatorial plane, the accretion flow is inclined to undergo a larger range 
of 6 to change Vr to 0, resulting in a larger inflow region. For the two SSD cases, as their flow 
structures are quite similar, larger values of Vr in the equatorial plane give the gas pressure dom- 
inated solution set ai larger inflow regions than a2. For the advection-dominated cases, not only 
the equatorial value of is smaller in the slim disk case than that in the ADAF case, but its radial 
direction pressure gradient also increases as 6 decreases, giving a larger outward acceleration to Vr, 
thus the gas pressure dominated solutions have larger inflow regions than the radiation pressure 
dominated solutions 04. 

The solution set 04, which corresponds to the slim disk model, has a large outflow size over a 
broad range of a, displaying a large outflow area. The solutions corresponding to the SSD model, 
ai and a2, have very small outflow regions, and their outflows are mainly in the form of wind 
blowing out of the disk surface instead of escaping in the outflow region. The solution set a^, 
which corresponds to the ADAF model, has a moderate outflow region with small a. However as 
a increases, this outflow region quickly disappears and the outflow takes the form of wind blowing 
out of the upper surface. The velocity field displayed in Figure 5 is an example for or = 0.1. 
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a is connected with both the viscous heating and the angular momentum transport of the ac- 
cretion flow. Larger a will increase the viscous heating for a fixed surface density. But at the same 
time it also increases the angular momentum transfer rate, thus increasing and decreasing surface 
density, which decreases the total viscous heating rate. Therefore the dependence of the viscous 
heating rate on a, is somewhat complex. Viscous heating represents how much gravitational en- 
ergy is converted to internal energy. But some of it is lost due to radiation and eventually there is 
only an / fraction of the viscous heating which is converted to internal energy. The increase in the 
internal energy will raise the temperature and subsequently the pressure, thus increasing the size of 
the accretion flow. However, for the SSD cases shown in ai & a2, f is very small, thus the energy 
converted to internal energy is negligible and the disk structure is dominated dynamically by the 
gravity and the rotation of the material, rather than the pressure gradient. As a result, the sizes of 
accretion flows change very little with diff'erent a for ai & 02- 

On the other hand, for the advection dominated solutions & 04, a large fraction of the 
gravitational energy is eventually stored as internal energy, so that the pressure gradient plays a 
more important role in determining the flow structure as mentioned in §3.1. Therefore the sizes of 
accretion flows are closely connected with a. In Figure 8 we can see that for 03 & 04, when a is 
very small, the viscous process is not quite effective, therefore both the inflow region size and the 
outflow region size are small (it's easy to understand if we consider the situation or = 0, in which 
case no accretion flow will be formed). From this small value, along with the increasing of a, 
the viscous process starts to take efi'ect and the whole inflow/outflow region size starts to increase, 
until it reaches a maximum size at or ~ 0.05 for line a'^ and a ~ 0.12 for line a'^(the exact value 
of a for maximum size also depends on other parameters). From then on, the larger values of 
starts to dominate the total efi'ect and the corresponding accretion flow size starts to decrease. This 
property can also be seen in the top left panel of Figure 17, which displays the ratio of the outflow 
rate to the inflow rate. 

It should be noted that, in Figure 8, the size change based on a is different between the whole 
inflow/outflow region and the pure inflow region. The size of the whole inflow/outflow region 
decreases after a peak at a certain a while the size of the pure inflow region keeps increasing. 
For solution sets ai, and a^, this difference is enough to generate a critical respectively, at 
which the outflow region totally disappears. So for a > a^, the accretion flow that we can resolve 
with self-similar assumptions is composed of pure inflow and we didn't find a surface at which 
VriOo) = 0. That's why we say that the assumption used in XW05 may not be true for large a. 
Kluzniak and Kita (2000) also found this critical a which was calculated to be V15/32 » 0.685 
by them. In our work, however, depends on other parameters and can change significantly, as 
shown in Figure 8 {~ 0.96 for a\, ~ 0.24 for and ~ 0.74 for 04). For 02, this ac is above 1 and 
it's not displayed. 



-15- 



Figures 9 displays the Vr distribution along the ^-direction for solutions with different a. Fig- 
ures (a), (b), (c) and (d) correspond to the four sets of solutions in Figure 8 respectively, and the 
solid, dashed and dotted lines correspond to different a as indicated in the legends. The range of 
the axes are adjusted to show more details in comparison. We can see that in the solutions cor- 
responding to the SSD model, i.e. Figures 9(a) & 9(b), Vr is almost proportional to a. That is 
because when alpha increases, so does the angular momentum transfer rate. For the SSD case, the 
accretion flow structure is mostly determined by gravity and rotation as discussed in §3.1, and the 
rotation is quasi-Keplerian, so the increase in angular momentum transfer rate increases almost 
proportionally. A special case in the equatorial plane is shown by Equation (31), in which when 
V0 is Keplerian, the value of Vr is proportional to a. For the ADAF and the slim disk cases, be- 
cause pressure gradient plays an important role in the accretion flow structure, Vr is generally not 
proportional to a, but still positively correlated with a as shown in Figure 9(c) & 9(d). We can 
also see that some lines in Figure 9 don't contain positive v^. That's because the a values for these 
solutions are greater than the critical value ac mentioned above. 

2) Solution dependence on n 

The parameter n describes how the density changes along the radius in our self-similar as- 
sumptions. As discussed in §3.1, when n = 3/2, the entire flow structure can be described by the 

same set of self-similar solutions, but doesn't contain outflow; when n < 3/2, the solution con- 
tains outflow but the self-similar assumptions cannot be applied to the region near the vertical axis; 
solutions with n > 3/2 are unlikely to happen in real cases. 

The dependence of 9o and ^ on n is shown in Figure 11. Figures 12 & 13 display the Vr and 
V0 distributions along the ^?-direction for solutions with different n. We can see that solutions all 
achieve a maximum size at n = 3 /2 for both the whole inflow/outflow region and the inflow only 
region. This is because when n = 3/2 the entire space can be described by the same set of self- 
similar equations as discussed in §3.1. In this case, = and mass is conserved on the streamlines 
pointing exactly at the central accretor, thus no outflow is required. Before reaching the maximum 
value atn = 3/2, the flow size is positively correlated with n as a general trend. As n decreases, 
according to Equation (22) M^ff decreases faster inwards, which means a larger mass loss rate into 
the area beyond 0^, resulting in the increase in the size of the wind region and the decrease in the 
size of the whole inflow/outflow region in Figure 11, and the increase in vg in Figure 13. 

In Figure 1 1 line b'^ has 2 peaks. According to Equation (25), the radial pressure gradient is 
proportional to (1 +n)p/p, so when p/p doesn't change much, the radial pressure gradient increases 
with n. In advection-dominated solutions & &4, the radial pressure gradient is more important 
than centrifugal force in accelerating v^, so for larger n which is inclined to have larger radial 
pressure gradient, one would expect to see faster increase in v,- and stronger outflow in the outflow 
region. These features of Vr can be seen in Figure 12. The radial pressure gradient is dominant 
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in the slim disk case as mentioned in §3.1, so the acceleration of Vr increases significantly with 
increasing n in the slim disk case, which can be seen in Figure 12(d) with a much steeper profile of 
Vr when n = 1.3. This increased outflow velocity in the outflow region can cause a steeper density 
drop along 9 direction towards the upper boundary and therefore make the outflow region become 
smaller, creating a valley in line b'^ around n = 1 .4. 

The top right panel in Figure 17 shows how the ratio of the outflow rate to the inflow rate 
changes depending on n. We can see that when n becomes small enough, the ratio changes very 
little with n. In this range of n, the size of different regions also doesn't change much in Figure 10, 
and the flow structure is similar to each other. 

For yequ = 5/3, we cannot get a reasonable set of equatorial values for physical properties 
when n is larger than 3/2. For ygqu = 4/3, i.e. solutions b2 and Z?4, we can still get reasonable values 
of equatorial physical properties when n is larger than 3/2, so that we can calculate solutions with 
n > 3/2 in this case. Solutions with n > 3/2 have an effective accretion rate Mes that increases 
towards the central accretor, which implies that there must be matter injection into the accretion 
flow from high latitudes. This kind of matter injection is not likely to happen in real cases, so we 
won't pay more attention to the solutions with n > 3/2 here. Solutions with n = 3/2 and n <3/2 
both exist and have quite different structure. However, as in NY95, solutions with n = 3/2 are 
only for ADAFs, while solutions with n < 3/2 in our work is applicable to a much wider range 
of accretion disks. In this sense, we believe that outflows should be common in various accretion 
disks. 

3) Solution dependence on / 

The dependence of and ^ on / is shown in Figure 14. It has already been discussed that as 
the advective factor / increases, the size of the whole inflow/outflow region also increases, which 
is shown more clearly here in Figure 14. The gas pressure dominated solutions have larger inflow 
region in agreement with Figures 8 & 11, and the reason has been discussed in §3.2.1. 

As / increases, the size of the outflow region increases faster for radiation pressure dominated 
flows, displaying a larger outflow size for solutions corresponding to the slim disk model. The 
reason is that, as / increases, a larger fraction of viscous heating is converted to internal energy 
instead of being lost via radiation, raising the temperature of the accretion flow. The gas pressure 
is proportional to T while the radiation pressure is proportional to T^, so in the radiation pressure 
dominated solutions, p increases faster as / increases. The increased pressure inflates the outflow 
region and causes a stronger outflow. On the other hand, as / gets close to 1, eventually the 
outflow velocity at small inclination becomes so large, that it will instead reduce the size of the 
outflow region because the material lost in the outflow region is limited. As shown in the bottom 
left panel of Figure 17, at a certain / ~ 0.9, the fraction of material lost in the outflow region 
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reaches an upper limit, and with larger /, this fraction almost remains unchanged. In this range of 
/, the larger outflow velocity will instead shrink the size of the outflow region. This causes line 
C2 in Figure 14 to reach a peak of ~ 90° at / ~ 0.9 (corresponding to Figure 17(c)), and to start 
decreasing as / becomes larger. 

4) Solution dependence on yequ (and consequently p) 

The dependence of and on yequ is shown in Figure 15. The top axis displays the corre- 
sponding values of yS when 7 takes the value of 5/3. If the intrinsic heat capacity ratio is not 5/3 
(e.g. when dealing with proto-stellar accretion disks with unionized gas), the corresponding p in 
Figure 15 may not be correct, but one can still calculate the correct p with yequ and y according 
to Equation (14). Figure 16 displays the and vg distributions along the ^-direction for solutions 
with different yequ- The bottom right panel of Figure 17 shows how the ratio of the outflow rate to 
the inflow rate changes depending on yequ (and consequently P). 

In Figure 15, the size of the pure inflow region generally increases with yequ (and conse- 
quently P). As shown in Figure 10 & 16 and discussed in §3.2.1, solutions with larger yequ (and 
consequently larger gas pressure ratio P) have larger values of v,. in the equatorial plane. This eff"ect 
combined with the acceleration properties of v,., causes solutions with larger yequ (and consequently 
larger P) to have larger inflow region. For solutions with low advection, as the outflow region is 
very small, the size of the whole inflow/outflow region also displays similar profile with increasing 
yequ (and consequently p). That's why lines di, d[ and d2 all have positive slope in Figure 15. 

In Figure 17(d), there exists a transition point for the advection dominated solutions, and for 
solutions with yequ < 3/2 (corresponding to < 2/3), the ratio of the outflow rate to the inflow 
rate keeps being large. As discussed before, for the advection dominated cases, radiation pressure 
dominated solutions have larger pressure gradient, shown as steeper velocity profiles in Figure 16, 
which drives stronger outflow. As shown in Figure 17(d), the outflow to inflow ratio increases as 
yequ decreases from 5/3 (which indicates the decreasing of gas pressure ratio p from 1). However, 
the outflow rate cannot exceed the inflow rate and eventually the increase will cease, reaching 
a transition point at yequ ~ 3/2 (corresponding to /3 ~ 2/3). For yequ < 3/2 (corresponding to 
yS < 2/3), the outflow to inflow ratio remains almost unchanged. Further decrease in yequ will 
continue to increase the outflow velocity, while the fraction of the inflow lost in the outflow region 
is almost unchanged, so the size of the outflow region will shrink. This feature can be also seen in 
Figure 15, which shows that the size of the outflow region for advection dominated solutions first 
increases with decreasing yequ, reaching an upper limit at yequ ~ 3/2 (corresponding to p ~ 2/3), 
and decreases with decreasing yequ after that. It should be noted that, this transition value of yequ 
depends on other input parameters, and this value of ~ 3/2 corresponds to or = 0.1, n = 1.3 and 
/ = 1 . For the solutions with low advection, the ratio of the outflow rate to the inflow rate don't 
change much and are always very small as shown in Figure 17(d). 
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5) Bends in solutions 

We have shown the solution dependence on different input parameters. The solutions some- 
times show bends which we would like to discuss more specifically here. 

The bend about a is closely connected with the energy mechanism. As discussed in §3.2.1, 
the dependence of viscous heating on a is somewhat complex. It has a peak at a certain value of a 
(which depends on other input parameters), and decreases on both sides. For the advection dom- 
inated solutions fl3 & a4, most of the viscous heating is stored as internal energy in the accretion 
flow, which raises the temperature and subsequently the pressure, thus increasing the size of the 
accretion flow. On the other hand, the increase in a also increases the absolute value of Vr of the 
inflow near the equatorial plane, in which case more work needs to be done to drive an outflow, de- 
creasing the outflow region as a increases. So naturally the size of the whole inflow/outflow region 
also has a peak at some value of a. As a increases from a very small value, at first the increase 
in viscous heating dominates the total efi'ect, and the size of the accretion flow increases, until it 
reaches a maximum value; from then on, the increase in Vr starts to dominate the total eff'ect, which 
makes the outflow more difficult to produce and also reduces viscous heating when a becomes 
large enough, thus shrinking the size of the whole accretion flow. To show the energy profile for 
solutions with difl'erent a, here we calculate the 0-averaged Bernoulli function Be. The Bernoulli 
function represents the specific total energy in the accretion flow. A positive value of the Bernoulli 
function is only a necessary, not a sufficient, condition for outflow formation (Abramowicz et al. 
2000). However, as our solutions have already showed outflow structure in the velocity fields, a 
larger value of the Bernoulli function indicates that the accretion flow is more likely to contain 
stronger outflow. Locally the Bernoulli function can be written as 



1 , P 1 9 GM 



i(v^(e) + v^(e) + v5(e)) + ^^^-l 



, (34) 



in which W is the specific enthalpy, V is the velocity (all three components included) and O is the 
gravitational potential. So the 0-averaged Bernoulli function can be written as 

_ C" p{e) ■ Be ■ Inr^ sin Odr ■ {nj 1 80°) J0 



Be = 



p(e) ■ 2nr^ sin Odr ■ (n/ 1 80°)^^^ 



Tequ 



-1 



sin Ode 



f p(e)-smede 



(35) 



Figure 18 displays how the 0-averaged Bernoulli function Be changes with variant a for the advec- 
tion dominated cases. It can be seen that the positions of the peaks agree quite well with the bends 
of line a'g and line a'^ in Figure 8. 
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Figure 19 displays the velocity fields of solutions around the bend for variant a. The values 
of a are chosen so that the solutions on both sides of the peaks have similar upper boundaries 
(~ 10°). It can be seen that the velocities in the outflow region are more aligned to the 6 direction 
in the right panel of each row than those in the left panel. That's because larger a increases the 
angular momentum transfer rate, thus increases the absolute value of Vr near the equatorial plane 
where most of the material is accreted. The value of Vy is negative in the equatorial plane and 
changes gradually to positive values in the outflow region under the eff'ect of pressure gradient and 
centrifugal force. The absolute value of Vy in the equatorial plane in the right panel of each row is 
significantly larger than that in the left panel, so when pressure gradient doesn't differ much, the 
positive values of Vy in the outflow region will be smaller for larger a and the velocities will be 
more aligned to the 6 direction. The solutions corresponding to slim disks (the second row) show 
velocities more aligned to the radial direction in the outflow region than those corresponding to 
ADAFs (the first row), due to larger pressure gradient and thus larger values of Vy in the outflow 
region. 

The bends in solution dependence on n, f and yequ have different physical mechanism from 
that on a. It has been mentioned that, as n increases, / increases or ygqu decreases, the pressure 
gradient increases and drives stronger outflow. However, the outflow rate cannot exceed the inflow 
rate. At seen in Figure 17, the ratio of the outflow rate to the inflow rate for lines b^, and 
will eventually reach an upper limit at certain values of corresponding parameters. For solution 
sets C2 and d^, further increase in / or further decrease in yequ will continue to increase the radial 
outflow velocity, while the outflow rate in the outflow region is almost unchanged, so the size of the 
outflow region will shrink. For the solution dependence on n, the bend near n = 1.25 has similar 
physical reason as the bends for / and yequ- However, the size of the accretion flow has another 
peak at n = 3 /2, where the flow is always radial (with rotation) and no outflow is needed. So line 
Z74 in Figure 17 starts to drop as n gets close to 3/2 and line h\ in Figure 1 1 gets to another peak at 
n = 3/2. 

Figure 20 displays the velocity fields of solutions around the bend for variant n, f and yequ- 
It can be seen that, as n increases (the first row), / increases (the second row) or yequ decreases 
(the third row), the flow patterns display similar evolvement for these parameters: at first both the 
accretion flow size and the radial outflow velocity increase; after the corresponding parameter gets 
past the respective value at the transition point (the middle panels in all three rows of Figure 20), 
the radial outflow velocity continues to increase while the size of the accretion flow decreases. The 
velocities in the outflow regions become more aligned to the radial direction from left to right in 
all three rows, due to increased radial outflow velocities. 

The labels for solution sets and lines which appear in the text and figures are summarized in 
Table 1 along with their corresponding input parameters. 
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Table 1 

Labels and corresponding parameters 



Label for solution set 


Label for lines 


a 




n 


/ 


Tequ 


ai 


ai, a[ 


0.01 - 


1 


1.3 


0.01 


5/3 


ai 


ai, a'2 


0.01 - 


1 


1.3 


0.01 


4/3 


as 


as, flj 


0.01 - 


1 


1.3 


1 


5/3 


a4 


04, «4 


0.01- 


1 


1.3 


1 


4/3 


bi 


bu b\ 


0.1 




0-1.5 


0.01 


5/3 


bi 


bi, b'2 


0.1 




0-1.5 


0.01 


4/3 


b3 


bi, ^3 


0.1 




0-1.5 


1 


5/3 


bA 


bA,b\ 


0.1 




0-1.5 


1 


4/3 


C\ 


Cu c\ 


0.1 




1.3 


0.01 - 1 


5/3 


C2 


Cl, 4 


0.1 




1.3 


0.01 - 1 


4/3 


di 


du d[ 


0.1 




1.3 


0.01 


4/3 - 5/3 


di 


d2,d'2 


0.1 




1.3 


1 


4/3 - 5/3 



4. SUMMARY AND DISCUSSION 

With the self-similar assumptions along the radial direction and boundary conditions obtained 
from the reflection symmetry in the equatorial plane, we are able to solve the Navier-Stokes Equa- 
tions along the 6 direction explicitly, and the velocity field is obtained. The result shows that 
outflows are common in all kinds of accretion disk models, which is consistent with some numer- 
ical simulation works (e.g. Ohsuga et al.2005; Ohsuga & Mineshige 2007; Ohsuga et al. 2009). 
Generally the accretion flow consists of three different regions: an inflow region near the equato- 
rial plane which contains the largest portion of mass, an outflow region above the inflow region in 
which the matter starts escaping the central accretor in the r-direction, and a wind region which 
contains the material blowing out from the boundary of the outflow region. The structure of the 
inflow and the outflow region can be resolved in our solutions, while the wind region doesn't obey 
the self-similar assumptions and can not be resolved by our calculation. The inflow region and the 
wind region are essential to form a steady accretion flow structure. The outflow region is missing 
in the cases with large viscosity parameter a > a^. The critical value depends on other parame- 
ters and is ~ 0.24 for ADAFs, ~ 0.74 for slim disks and ~ 0.96 for gas pressure dominated SSDs 
when n = 1.3. 

To compare with the popular analytical models, we calculated solutions with parameters cor- 
responding to the SSD, the ADAF and the slim disk models. The solutions corresponding to the 
SSD model have relatively thin, quasi-Keplerian disks with very small outflow regions, and the 
outflows mostly take the form of wind from the disk surface. The solutions corresponding to the 
ADAF and slim disk models have thick, sub-Keplerian disks. The solution corresponding to the 
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ADAF model has a moderate outflow region which shrinks quickly with the increase of a, and 
the outflow mainly takes the form of wind. However as the upper boundary of the outflow region 
is close to the vertical axis, eventually the outflow may take the form of bipolar jets, though we 
cannot actually resolve the explicit structure around the vertical axis. The solutions corresponding 
to the slim disk model has a large outflow region in which most of the outflow takes place, and 
near the upper boundary of the outflow region, the material escapes with supersonic velocity. 

Our calculation displays very small outflow regions for input parameters corresponding to the 
SSD model, and for these solutions the inflow region has similar structure with the one-dimensional 
SSD model. However the accretion rate inside the disk may not be conserved along the radius due 
to mass lost via wind. For the ADAF and the slim disk cases, our result displays an outflow region 
too large to be neglected, and we think it's necessary to consider the explicit structure in the or 
z direction to get more realistic results than the traditional analytical models. It was proposed in 
several works (Gu & Lu, 2007; Jiao et al. 2009) that the slim disk model has an upper limit of 
accretion rates, above which outflows seem to be inevitable. Later in another paper (Jiao & Lu 
2009), one-dimensional steady transonic global solutions for slim disks with presumed amount of 
outflows were calculated, and in this case the proposed upper limit of accretion rates for slim disks 
can be exceeded. Another paper written by Takeuchi et al. (2009) used the outflow data from two- 
dimensional radiation-hydrodynamic (RHD) simulations to construct an one-dimensional steady 
solution, and the result shows that the emergent spectra do not sensitively depend on the amount of 
mass outflow. However, these models are still one-dimensional models and they only incorporated 
the mass loss effect of outflow into their models, so the results still remain debatable. We think that 
two-dimensional solutions considering the hydrodynamical processes and radiative transfer along 
both the radial and the direction are required to examine the topic, which we are planning to do 
in our future work. 

We also investigated the dependence of the accretion flow structure on diff'erent input param- 
eters respectively. / and y^qu (consequently yS) both have some transition values as discussed in 
§3.2.3 and §3.2.4, at which the outflow region reaches a maximum size. Generally when / ^ 0.9 
and y8 < 2/3 (for intrinsic y = 5/3), the accretion flow will have a large outflow region and strong 
outflow, as the aforementioned solutions corresponding to the slim disk model. However, in what 
conditions the transition happens depends on other input parameters and it won't happen in case of 
large a due to the existence of the critical value qtc. Solutions with n = 3/2 don't contain outflow 
and they correspond to the solutions in NY95, which is only applicable to ADAFs. Solutions with 
n < 3/2 have the aforementioned three-region (or two-region for a > ac) structure. Solutions with 
n > 3/2 are not likely to happen in real cases. 

Finally, we would like to mention some caveats in our work. The solutions obtained in this 
paper are steady solutions, and thus it is not guaranteed that they are stable solutions. These solu- 
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tions do not cover the whole space, so for the unresolved region near the axis, we know nothing 
about the detailed structure other than that it contains the material blowing out of the upper bound- 
ary in the form of wind. Our work is based on the simple ap prescription of viscosity, which 
is not bad according to recent MHD simulation works(Hirose et al. 2009; Ohsuga et al. 2009), 
especially when the flow is steady. However, it may be necessary to consider other components 
of the viscous stress tensor for solutions with large v, . We also used the advective factor / to 
simplify the energy equation, and assumed / to be a constant. Realistically / should vary with 
both 6 and r, which could be improved by considering the details of radiative transfer. We used 
Newtonian gravitational potential in our work, and it needs to be corrected when applied to study 
the structure close to the innermost stable circular orbit. In our calculation we didn't consider 
the eff"ects of convection (e.g. Stone et al. 1999; Igumenshchev & Abramowicz 2000; Yuan & 
Bu 2010) or magnetic field (e.g. Blandford & Payne 1982; Naso & Miller 2010), which may be 
important in studying the accretion flow and the generation of outflows. Finally, it is necessary to 
abandon the self-similar assumptions if we want to investigate the disk and outflow structure more 
quantitatively and precisely. These caveats will be improved in our future work. 
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his/her constructive suggestions. This work was supported by the NSFC grant (No. 11033001) 
and the 973 program (No. 2007CB8 15405). 

REFERENCES 

Abramowicz, M. A., Chen, X., Kato, S., Lasota, J.-R, & Regev, O. 1995, ApJ, 438, L37 
Abramowicz, M. A., Czemy, B., Lasota, J. R, & Szuszkiewicz, E. 1988, ApJ, 332, 646 
Abramowicz, M. A., Lasota, J.-R, & Igumenshchev, I. V. 2000, MNRAS, 314, 775 
Blandford, R. D., & Begelman, M. C. 1999, MNRAS, 303, LI 
Blandford, R. D., & Begelman, M. C. 2004, MNRAS, 349, 68 
Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 
Chartas, G., Brandt, W. N., & GaUagher, S. C. 2003, ApJ, 595, 85 
Dotan, C, & Shaviv, N. J. 2010, arXiv: 1004. 17971 
Fukue, J. 1989, PASJ, 41, 123 



-23- 



Fukue, J. 2004, PAS J, 56, 181 

Gu, W.-M., & Lu, J.-F. 2007, ApJ, 660, 541 

Hirose, S., Blaes, O., & Krolik, J. H. 2009, ApJ, 704, 781 

Igumenshchev, I. V., & Abramowicz, M. A. 2000, ApJS, 130, 463 

Jiao, C.-L., & Lu, J.-F. 2009, Chinese Physics Letters, 26, 049701 

Jiao, C.-L., Xue, L., Gu, W.-M., & Lu, J.-F 2009, ApJ, 693, 670 

Kato, S., Fukue, J., & Mineshige, S. 2008, Black-Hole Accretion Disks: Towards a New Paradigm 
(Kyoto: Kyoto Univ. Press) 

Kawabata, R., & Mineshige, S. 2009, PASJ, 61, 1135 

King, A. R., Pringle, J. E., & Livio, M. 2007, MNRAS, 376, 1740 



Kluzniak, W., & Kita, D. 2000, arXiv:astro-ph/0006266 



Kusunose, M. 1991, ApJ, 370, 505 

Loeb, A., Narayan, R., & Raymond, J. C. 2001, ApJ, 547, L151 
Marrone, D. P, Moran, J. M., Zhao, J.-H., & Rao, R. 2006, ApJ, 640, 308 
Narayan, R., & Yi, I. 1994, ApJ, 428, L13 
Misra, R., & Taam, R. E. 2001, ApJ, 553, 978 
Narayan, R., & Yi, I. 1995, ApJ, 444, 231 
Naso, L., & Miller, J. C. 2010, A&A, 521, A31 

Okuda, T., Teresi, V., Toscano, E., & Molteni, D. 2005, MNRAS, 357, 295 
Ohsuga, K., Mori, M., Nakamoto, T., & Mineshige, S. 2005, ApJ, 628, 368 
Ohsuga, K., & Mineshige, S. 2007, ApJ, 670, 1283 
Ohsuga, K., Mineshige, S., Mori, M., & Kato, Y. 2009, PASJ, 61, L7 

S^dowski, A., Abramowicz, M., Bursa, M., Kluzniak, W., Lasota, J. -P., & Rozahska, A. 2011, 
A&A, 527, A17 

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 



Stone, J. M., Pringle, J. E., & Begelman, M. C. 1999, MNRAS, 310, 1002 

Takahara, R, Rosner, R., & Kusunose, M. 1989, ApJ, 346, 122 

Takeuchi, S., Mineshige, S., & Ohsuga, K. 2009, PASJ, 61, 783 

Xie, F.-G., & Yuan, F. 2008, ApJ, 681, 499 

Xu, G., & Chen, X. 1997, ApJ, 489, L29 

Xue, L., & Wang, J. 2005, ApJ, 623, 372 

Yuan, R, & Bu, D.-R 2010, MNRAS, 408, 1051 



This preprint was prepared with the AAS U^^i macros v5.2. 



-25- 




Fig. 1. — The solution corresponding to the gas pressure dominated region of the SSD model. 
Here a = 0.1, n = 1.3, / = 0.01 and yequ = 5/3, which correspond to gas pressure dominated 
monatomic ideal gas. 




Fig. 2. — The solution corresponding to the radiation pressure dominated region of the SSD model. 
Here a = 0.1, n = 1.3, f = 0.01 and yequ = 4/3, which correspond to radiation pressure dominated 
monatomic ideal gas. 
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Fig. 3. — The solution corresponding to the ADAF model. Here a = 0.1, n = 1.3, / = 1 and 
7equ = 5/3, which correspond to gas pressure dominated monatomic ideal gas. 
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Fig. 4. — The solution corresponding to the sUm disk model. Here a = 0.1, n = 1.3, / = 1 and 
7equ = 4/3, which correspond to radiation pressure dominated monatomic ideal gas. 
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Fig. 5. — The velocity fields of the four typical solutions. Figures (a), (b), (c) and (d) correspond 
to input parameters of Figures 1-4 respectively. The lengths of arrows indicate the absolute val- 
ues of the vector v',-(d)+Vf){9), which are scaled logarithmically. The solid lines correspond to the 
inclination 6^,, while the dashed lines correspond to the inclination ^o- 
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Fig. 6. — The radial components of the pressure gradient, the centrifugal force and the gravitational 
force for the four typical solutions. Figures (a), (b), (c) and (d) correspond to input parameters of 
Figures 1-4 respectively. The dotted lines correspond to the gravitational force which is scaled 
as 1; the dashed lines correspond to the radial component of the centrifugal force; the dot-dashed 
lines correspond to the radial component of the pressure gradient; and the solid lines correspond to 
the sum of the radial components of the centrifugal force and the pressure gradient, which drives 
the outflow. 
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Fig. 7. — The solution corresponding to et = 0.1 in NY95. Here a = 0.1, n = 1.5, / = 1 and 
7equ = 1.6061. 
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Fig. 8. — Solution dependence on a. Solid lines ai, ai, 03 and 04 represent the upper boundaries 
of the inflow region, and dotted lines a\, a'^, a'^ and a\ represent the boundaries of the whole 
inflow/outflow region. All the lines correspond ion = 1.3. Lines a\ and a'j correspond to / = 0.01 
and 7equ = 5/3; lines ai and a\ correspond to / = 0.01 and yequ = 4/3; lines as and 03 correspond 
to / = 1 and 7equ = 5/3; lines 04 and a\ correspond to / = 1 and yequ = 4/3. 
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Fig. 9. — The distribution of Vr(0) along the ^-direction for solutions with different a. Figure (a) 
corresponds to n = 1.3, f = 0.01 and yequ = 5/3; Figure (b) corresponds to n = 1.3, f = 0.01 and 
7equ =4/3; Figure (c) corresponds to n = 1.3, f = 1 and yequ = 5/3; Figure (d) corresponds to 
n = 1.3, / = 1 and yequ = 4/3. The solid, dashed and dotted lines correspond to or = 0.1, 0.5 and 
1, respectively. 
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Fig. 10. — The equatorial values of vXd) corresponding to different yequ- AH the solutions cor- 
respond to a = 0.1 and n = 1.3. The solid line corresponds to / = 1, while the dashed line 
corresponds to / = 0.01. on the top axis corresponds to yequ on the bottom axis and is calculated 
based on y = 5/3. 
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Fig. 11. — Solution dependence on n. Solid lines bi, b2, and b4 represent the upper boundary 
of the inflow region, and dotted lines b[, b'2, b'^ and b'^ represents the boundaries of the whole 
inflow/outflow region. All the lines correspond to a = 0.1. Lines bi and b\ correspond to / = 0.01 
and 7equ = 5/3; lines &2 and b'^ correspond to / = 0.01 and yequ =4/3; lines b^, and b'^ correspond 
to / = 1 and 7equ = 5/3; lines b^. and b'^ correspond to / = 1 and yequ = 4/3. 
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Fig. 12. — The distribution of Vr(0) along the 0-direction for solutions with different n. Figure (a) 
corresponds to a = 0.1, / = 0.01 and 7equ = 5/3; Figure (b) corresponds to a = 0.1, / = 0.01 and 
7equ = 4/3; Figure (c) corresponds to or = 0.1, / = 1 and yequ = 5/3; Figure (d) corresponds to 
a = 0.1, / = 1 and 7equ = 4/3. The solid, dashed and dotted lines correspond to n = 1.3, 1 and 0, 
respectively. 
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Fig. 13. — The distribution of vo(6) along the 0-direction for solutions with different n. Figure (a) 
corresponds to o- = 0.1, / = 0.01 and ygqu = 5/3; Figure (b) corresponds to o- = 0.1, / = 0.01 and 
Tequ = 4/3; Figure (c) corresponds to a = 0.1, / = 1 and ygqu = 5/3; Figure (d) corresponds to 
a = 0.1, / = 1 and ygqu = 4/3. The solid, dashed and dotted lines correspond to n = 1.3, 1 and 0, 
respectively. 
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Fig. 14. — Solution dependence on /. Solid lines ci & C2 represent the upper boundary of the 
inflow region, and dotted lines c\ & c'2 represents the outer boundary of the whole accretion flow. 
All the lines correspond to o' = 0.1 & n = 1.3. Lines Ci and c[ correspond to yequ = 5/3; lines C2 
and c'2 correspond to 7equ =4/3. 
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Fig. 15. — Solution dependence on yequ- Solid lines di & J2 represent the upper boundary of the 
inflow region, and dotted lines d[ & d'2 represent the upper boundary of the whole inflow/outflow 
region. All the lines correspond to or = 0.1 & n = 1.3. Lines di and d[ correspond to / = 0.01; 
lines d2 and d'2 correspond to / = 1. jS on the top axis corresponds to yequ on the bottom axis and is 
calculated based on 7 = 5/3. 
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Fig. 16. — The Vr and vg distributions along the ^-direction for solutions with different yequ- Figures 
(a) & (c) correspond to solutions with n = 1.3, a = 0.1 and / = 0.01, while Figures (b) & (d) 
correspond to solutions with n = 1.3, a = 0.1 and / = 1. The solid, dashed and dotted lines 
correspond to yequ = 4/3, 3/2 and 5/3, respectively. 
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Fig. 17. — The ratio of the outflow rate to the inflow rate. The vertical axes show the absolute 
values of the ratio, and the horizontal axes correspond to different input parameters. The lines in 
the top left panel correspond to the lines in Figure 8. The lines in the top right panel correspond to 
the lines in Figure 11. The lines in the bottom left panel correspond to the lines in Figure 14. The 
lines in the bottom right panel correspond to the lines in Figure 15. Certain parts of some lines are 
missing due to that the solutions there don't contain an outflow region. 
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Fig. 18. — The ^-averaged Bernoulli function for variant a. Line as corresponds to solutions with 
/ = 1, 7equ = 5/3 and n = 1.3. Line ^4 corresponds to solutions with / = 1, 7equ = 4/3 and 
n = 1.3. 
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Fig. 19. — The velocity field patterns around the bends for variant a. The first row of figures 
corresponds to the bend of line a'^ in Figure 8, with a = 0.043, 0.05 and 0.12 from left to right 
respectively, and other parameters / = 1, 7equ = 5/3 and n = 1.3. The second row of figures 
corresponds to the bend of line a'^ in Figure 8, with a = 0.106, 0.12 and 0.165 from left to right 
respectively, and other parameters / = 1, 7equ = 4/3 and n = 1.3. The lengths of arrows indicate 
the absolute values of the vector v^(0)+v^(0), which are scaled logarithmically. The solid lines 
correspond to the inclination ^, while the dashed lines correspond to the inclination ^o- 
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Fig. 20. — The velocity field patterns around the bends for variant n, f and yequ- The first row 
of figures corresponds to the bend of line b'^ in Figure 11, with n = 1.1, 1.25 and 1.4 from left 
to right respectively, and other parameters a = 0.1, / = 1 and 7equ = 4/3. The second row of 
figures corresponds to the bend of line in Figure 14, with / = 0.8, 0.9 and 1 from left to right 
respectively, and other parameters a = 0.1, yequ = 4/3 and n = 1.3. The third row of figures 
corresponds to the bend of line d'2 in Figure 15, with yequ = 47/30, 3/2 and 43/30 from left to right 
respectively, and other parameters or = 0.1, / = 1 and n = 1.3. The lengths of arrows indicate 
the absolute values of the vector vv(^)+ve(6'), which are scaled logarithmically. The solid lines 
correspond to the inclination 9^, while the dashed lines correspond to the inclination ^o- 



